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We theoretically study bound states generated by magnetic impurities within conventional s-wave 
superconductors, both analytically and numerically. In determining the effect of the hybridization of 
two such bound states on the energy spectrum as a function of magnetic exchange coupling, relative 
angle of magnetization, and distance between impurities, we find that quantum phase transitions 
can be modulated by each of these parameters. Accompanying such transitions, there is a change 
in the preferred spin configuration of the impurities. Although the interaction between the impu¬ 
rity spins is overwhelmingly dominated by the quasiparticle contribution, the ground state of the 
system is determined by the bound state energies. Self-consistently calculating the superconducting 
order parameter, we find a discontinuity when the system undergoes a quantum phase transition as 
indicated by the bound state energies. 
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Introduction .—In a conventional s-wave superconduc¬ 
tor, quasiparticle excitation energies are separated from 
the chemical potential due to the formation of the super¬ 
conducting gap. When magnetic impurities are present, 
the exchange interaction can induce a bound state within 
the gap known as a Yu-Shiba-Rusinov (YSR) state, ^ 
which has been studied in detail both experimentally and 
theoretically.Recently, these states have attracted 
much attention in the context of magnetic impurity 
chains in which, when sufficiently close together, individ¬ 
ual YSR states can hybridize with adjacent bound states 
forming a band within the superconducting gap that can 
host Majorana fermions at its ends.^®“^® 

Two magnetic impurities interacting via quasiparti¬ 
cles are well described by the Ruderman-Kittel-Kasuya- 
Yosida (RKKY) interaction^^^^^ when the exchange in¬ 
teraction between the impurity and quasiparticles is 
much smaller than the Fermi energy. This results in a 
noncollinear orientation between the impurities in three- 
dimensional superconductors.Although for many 
parameters the contribution to the inter-impurity ex¬ 
change mediated by the overlap of YSR states is much 
smaller than that of the quasiparticles,it has been 
shown that resonant YSR bound states can dominate 
the exchange interaction and induce an antiferromagnetic 
alignment of the impurities.However, for the experi¬ 
mentally relevant limit^^ when the exchange interaction 
is of the order the Fermi energy, a theoretical under¬ 
standing of the interaction between magnetic impurities 
including (1) the quasiparticle contribution and (2) a self- 
consistent local reduction of the gap is missing from the 
literature. 

In this Letter, we determine the interaction between 
two magnetic impurities for arbitrary angles between 
their spins wherein the strength of the exchange interac¬ 
tion is unrestricted and, in general, unequal at the sites 
of the impurities. First, by analytically calculating the 
bound state energy spectrum, we find that a quantum 
phase transition (QPT)^“®d4 tuned by changing 



FIG. 1. Our setup of two magnetic impurities at ri and Y 2 
in an s-wave superconductor with classical spins Si and S 2 , 
respectively, oriented at a relative angle 0. As a result of the 
magnetic exchange couplings, Ji and J 2 , YSR bound states 
form within the bulk gap Aq. When the distance between 
the impurities, r, is larger than the coherence length of the 
superconductor the energies are Ei and E 2 but get changed 
to ei and 62 as r decreases and the bound states hybridize 
with each other. 


the distance between and relative magnetic orientation 
of the impurities. We, numerically, include the bulk 
contribution to the exchange interaction which quanti¬ 
tatively dominates over the YSR contribution for many 
parameters.Further, carrying out self-consistently 
calculations, we find a discontinuity in the superconduct¬ 
ing order parameter when the system undergoes such a 
QPT as indicated by the bound state energies. This, in 
turn, gives rise to magnetic metastable states, in addi¬ 
tion to the lowest energy magnetic configuration, for a 
sufficiently large exchange interaction. 

Model. We consider two magnetic impurities embed¬ 
ded in a bulk s-wave superconductor, see Fig. 1. The 
quasiparticles interact with the impurity spin through 
the exchange interaction that produces a local effective 
magnetic field. The corresponding Bogoliubov-de Gennes 
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Hamiltonian density is given by 

H = + A(r)Ta; - ^ JiSi ■ a-Sivi - r), (1) 

i=l,2 


where is the dispersion of the quasiparticles with mo¬ 
mentum p in the normal metal phase and A(r) is the 
local superconducting pairing strength. The Pauli ma¬ 
trices T (cr) act in Nambu (spin) space. The exchange 
coupling strength Ji of the spin impurity at can be pos¬ 
itive or negative corresponding to ferro- or antiferromag¬ 
netic interactions with quasiparticles, respectively. Here, 
we focus on Ji > 0 without loss of generality. We assume 
that Si are the classical spin vectors of the impurity at 
Ti, and 9 is the angle between them. The magnitude of 
the spins. Si = |Si|, are much larger than h so that quan¬ 
tum mechanical spin fluctuations, e.g. the Kondo effect, 
are negligible. In the following analytics we assume that 
A(r) = Aq is spatially uniform and neglect its suppres¬ 
sion due to the impurities,®’^® which we account for self- 
consistently in the numerics following earlier work.^“®’^^ 
To determine the energy of the bound states, € 1^2 and 
—€ 1,2 (particle-hole symmetry), we apply a straightfor¬ 
ward calculation along the lines of Ref. 21 and obtain 
a coupled set of secular equations for the BdG four- 
component spinors 'ip(r) at ri and r 2 . 


V'(i'i) = -^iSi • crV’(ri) -k f2S2 • crip(r2 ), 

V'(r2) = J2S2 • cr^(r2) -k f iSi • crip(ri ), 
where = Si/Si and 

f _ ai{e + TxAq) 

y, f{€ + TxAo)smkFr \ 


( 2 ) 

( 3 ) 

( 4 ) 


for f = 1, 2. Here, ai = VQTrJiSi, where vq is the density 
of states evaluated at the Fermi energy, r = |ri — r 2 | is 
the distance between the impurities, kp (vf) is the Fermi 
wave vector (Fermi velocity) and = vp/^ Aq — e^. 
When the distance between impurities is much greater 
than the superconducting coherence length, r ^ ^O) 
the impurities effectively decouple, F^ —)■ 0, and one 
finds that Eq. (2) furnishes solutions that are non¬ 
overlapping YSR bound states at ri and r 2 with ener¬ 
gies ±Ei = ±Ao(l — Q;f)/(1 -I- In this limit, for 

sufficiently large exchange interaction, Ji > l/v^pSi, the 
bound state energy goes below the chemical potential and 
the system undergoes a QPT wherein the parity of the 
ground changes. 

In order to determine the energies of the hybridized 
bound states analytically from Eq. (2), we focus on dis¬ 
tances between impurities much smaller than the coher¬ 
ence length, r <C ^ 0 , so that « 1 and the hybridiza¬ 

tion is determined to leading order by l/kpr. Formally 
diagonalizing the Hamiltonian and using a variational 


wave function as an ansatz for the ground state,®® the 
total energy of the system £gr is given by®’®® 

£A(^) = -\T.\^n{e)l (5) 

n 

where n, in general, runs over all solutions to Eq. (1); 
in the following analytics we only sum the bound state 
energies and determine £{9) = —(|ei(0)| -I- |e2(d)|)/2. 

Weak hybridization. For the moment, we consider the 
case of weak hybridization {kpr ^ Ei /Aq) for YSR 
states sufficiently far away from the chemical potential, 
so that the occupation of the bound states, and thus 
the ground state, is fixed by ai. That is, when ai < 1 
(ai > 1) the energy is above (below) the chemical po¬ 
tential. Calculating the full analytic solution and then 
expanding to second order in l/kpr, which is valid when 
|1 — ai\kpr 1 and \ai — a 2 \kpr 1,®® the spectrum 
has two solutions of the form 


e„( 6 ») K. En + Ao(H„ -k cos 9)/(kpr)"^, ( 6 ) 

where the coefficients and are functions of ai, a 2 , 
and kpr.^^ The bound state energy is extremized when 
either 0 = 0 or tt, i.e. the groundstate of impurities is 
collinear. When £162 > 0, £{'k) is always smaller than 
£( 0 )®®’®® and therefore the ground state is antiferromag¬ 
netic. When £162 < 0, £( 11 ) > £(0) and a ferromagnetic 
orientation is favored. 

Strong hybridization: identical impurities. Although 
strong hybridization between impurities cannot be ad¬ 
dressed perturbatively, in the symmetric case of equal ex¬ 
change coupling, i.e. ai = q ;2 = a, Eq. (2) can be solved 
directly. Because the analytic solution for arbitrary 9 is 
too involved, we focus here on collinear alignments. In 
the ferromagnetic configuration, the bound state energy 
levels are given by 

e± = eu2(0) = -Ao(a ± 6)/\/(a ± by + c® , (7) 


where a, b, and c, which depend on kp, r and a are dis¬ 
cussed in Ref. 38. The initially twofold degenerate energy 
levels of the bound states are both split due to hybridiza¬ 
tion and shifted due to the effective Zeeman splitting at 
both ri and r 2 . 

In the antiferromagnetic configuration, the energy level 
stays twofold degenerate® and is given by®® 


= ei,2(7r) = Aq 


I (1 — + 2{a/kFr)‘^ + d 

(1 + a®)® -k 2{a/kprY cos2kpr + d 


( 8 ) 


The difference in YSR bound state energy between 
the two collinear configurations, S£ = £(0) — f (tt) = 
+ I^fI ~ 2|£^|)/2, as a function of kpr is shown 
in Fig. 2. When a = 0.5 [Fig. 2 (upper panel)], all the 
electron-like energies in either configuration are greater 
than zero, in the displayed range, kpr > 1. Furthermore, 
S£ > 0 and therefore the exchange interaction between 
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FIG. 2. Energy of the YSR bound states for the identical mag¬ 
netic impurities oriented ferromagnetically (solid and dashed) 
and antiferromagnetically (dotted) as well as the energy dif¬ 
ference 5£ (thick solid lines) as a function of the distance 
r between impurities. When a = 0.5 (top panel), the sys¬ 
tem remains antiferromagnetic {S£ > 0), while for a = 0.9 
(lower panel) the magnetic configuration oscillates between 
being antiferromagnetic and ferromagnetic. For convenience, 
these two configurations are separated by the vertical dotted 
lines. 


impurities is antiferromagnetic, in agreement with the 
weak coupling limit. If the impurity levels are close to 
the chemical potential, e.g. a = 0.9 [Fig. 2 (lower panel)], 
the ground state of the system depends on the distance 
between the impurities. When r is sufficiently large, so 
that the condition for weak hybridization is met, the pre¬ 
ferred ordering is antiferromagnetic. When kpr 8, Cp 
goes below the chemical potential. Near this value of 
kpr, S£ becomes negative and therefore the preferred 
magnetic ground state is ferromagnetic rather than an¬ 
tiferromagnetic. As the distance between the impuri¬ 
ties decreases further, the bound state energies oscillate 
about the chemical potential as a function of r, thereby 
changing the YSR ground state. As a result, S£ also 
oscillates around zero implying a change between ferro¬ 
magnetic and antiferromagnetic configurations. 

Angle controlled quantum phase transition. As seen 
in the previous section, for some values of r (|Ai|/Ao 
kpr), the bound state energies are on opposite sides of 
the chemical potential in the ferromagnetic configuration 
due to hybridization, while in the antiferromagnetic the 
energies are always degenerate. Therefore, quite remark¬ 
ably, one may drive a QPT by changing the relative angle 
of the impurities. As shown in Fig. 3, one of YSR bound 
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FIG. 3. The energy of the bound states (ei, € 2 ) and the 
total YSR state energy (f (0)) as a function of relative angle 
9 for kpr = 1, ai = 0.5, and a 2 = 1. The change of quantum 
ground state at 0 « ±7^12 is indicated by vertical dotted lines. 

states passes through the chemical potential at 9 « 7r/2, 
signaling a QPT. The energy of YSR states is a minimum 
for the antiferromagnetic configuration, 9 = tt. Decreas¬ 
ing the angle between the impurities increases the energy 
until a critical point {9 « ±7r/2) when the ferromagnetic 
configuration becomes a minimum. Therefore, while the 
parameters chosen favor an antiferromagnetic configura¬ 
tion as an absolute ground state, they additionally sup¬ 
port a metastable ferromagnetic conhguration.'^® 

Bulk contribution. To address the contributions com¬ 
ing from the bulk, we follow earlier work^“®’^^ and 
study numerically a two-dimensional system with two 
magnetic impurities determining self-consistently the 
renormalization of the gap which cannot be addressed 
analytically. We use the tight-binding Hamiltonian 

H = -t Y + ^(A.CiiGi + H.c.) 

(T—±1 i 

+EE ([/X — 4:t + ((5ii + 5i2)JiCF cos9j\c\^Ciu 

i ( 7 — ±1 

T {dii bi 2 ^Ji sin ; (9) 

where Cia is the annihilation operator acting on an elec¬ 
tron with spin tr at lattice site i, and the first sum 
runs over neighboring sites i and i' located in a two- 
dimensional square lattice of size x Ny with lattice 
constant a. The chemical potential ^ is taken from the 
bottom of the energy band, and the local order parameter 
Ai is determined self-consistently in an iterative fash¬ 
ion for fixed values of the exchange coupling Ji at site 
i starting from the uniform superconducting order pa¬ 
rameter Aq. To compare to the analytics, we consider 
two impurities located at i = 1 and i = 2 (which are 
not necessarily adjacent) with equal exchange coupling, 
J = Ji = J 2 , and fixing the difference in magnetic orien¬ 
tation to be 9, mirroring the schematics of Fig. 1. After 
numerically diagonalizing Eq. (9), we find two types of 
energies in the spectrum: the energy of two YSR bound 
states £{9) considered before analytically and the total 
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FIG. 4. The difference in the energy 5£gr between ferro- and 
antiferromagnetic configurations of the system consisting of 
two identical impurities of coupling strength J as a function 
of the distance between impurities, r/a for (a) J/t = 1 and 
(b) J/t = 2.5 found self-consistently (red dots) and not self- 
consistently (blue dots). Insert: enlarged area of (a) for large 
distances. The difference in the energy 5S between ferro- 
and antiferromagnetic configurations including only the YSR 
bound state (green dots) is found self-consistently. The pa¬ 
rameters used are x Ny = 33x25, [ijt = 1, and Ao/t = 0.1. 


bulk energy £gr{0) obtained by summing all the energies 
below the chemical potential, see Eq. (5). 

First, we consider the difference between the ground 
state energies in the collinear magnetic configurations, 
6£ and d£gr as a function of distance r, see Fig. 4. The 
YSR bound state contribution S£ is positive for nearly all 
values of r, when Jjt = 1 [Fig. 4(a)], i.e. antiferromag¬ 
netic configuration is preferred. Whereas for J/t = 2.5 
[Fig. 4(b)], 6£ oscillates between positive and negative 
values. Both results agree with the analytics. 

Second, we aim to address the effect of gap renormal¬ 
ization and plot d£gr without self-consistent renormaliza¬ 
tion assuming = Aq (see Fig. 4). Interestingly, 5£gr 
is changed only slightly for all values of J, keeping the 
energies at the same order of magnitude. Upon including 
renormalization of the gap, 5£gr is increased drastically 
and the magnetic orientation becomes very sensitive to 
the distance between the impurities. This emphasizes the 
importance of a self-consistent renormalization of the gap 
when calculating the energies of such a system especially 
close to the phase transition. 

Third, we determine the angular dependence of the 
total energy and YSR bound state energy. For the fixed 
distance between the impurities (see Fig. 5), we observe 
that away from the phase transition, £gr and £ changes 
monotonically for 9 e [0,7r], and, thus, the ground state 


is either ferromagnetic or antiferromagnetic which is con¬ 
sistent with analytical results. In contrast to that, close 
to the phase transition when the bound state energies do 
cross the chemical potential as a function of 6, the depen¬ 
dence is non-monotonic [see Fig. 5(a)], and, in addition 
to the ferromagnetic (antiferromagnetic) ground state, 
there is a metastable antiferromagnetic (ferromagnetic) 
state. We also note that self-consistent solution demon¬ 
strates a jump in energy as one of YSR states crosses zero 
energy. Thus, we again find the qualitative agreement 
with analytical calculations predicting metastable states 
by analyzing only YSR bound states. However, we em¬ 
phasize that it is the QPT that results in the metastable 
state in Fig. 5(a), while the interaction is dominated by 
the bulk (not bound) state contribution to the energy. 


a) 
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FIG. 5. The total energy of the ground state for two identical 
impurities (J = Ji = J 2 ) as a function of the angle 9 between 
magnetic moments found numerically for (a) J/t = 2.17 (fer¬ 
romagnetic ordering) and (b) J/t = 2 (antiferromagnetic or¬ 
dering) at the distance r/a = 6. Other parameters are the 
same as in Fig. 4. 


Conclusions .—We have studied how the orientation of 
two spin impurity coupled via overlap of the YSR bound 
states induced by them depends on the distance between 
impurities and the strength of the exchange interaction. 
We have also demonstrated that a QPT can be controlled 
by changing relative magnetic orientation. Generally, the 
bulk contribution to the total ground state energy dom¬ 
inates over the bound state contribution, especially if 
the superconducting order parameter is determined self- 
consistently. The proposed effects could be measured 
with STM^® or NV-center^^ds techniques. 
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Appendix A: Variational Wave Function 


We extending the variational wave function for one impurity® to two impurities. For sufficiently weak coupling, in 
both the exchange interaction (ai, a 2 ^ 1) and the bound state hybridization {kFr ^ 1 ), the ground state is given by 
the BCS-like wave function |5'o) ~ nn>o('*^" + where ipn furnish a basis for the BdG Hamiltonian in the 

presence of the impurities for a given magnetic alignment and and are the Bogoliubuv coherence factors. The 
quasiparticle operators 7 „ are defined as 71 = — 7 !]^ = ui'tp-i+vi'ipl, 7 J = uiip\—viip-i, and similarly for 

n > 1, so that 7n|'I'o) = 0 for all n. Let n = 1 correspond to the lower energy bound state and n = 2 to the higher one 
while —n corresponds to a state with reversed spin. When the lower energy bound state is occupied, the wavefunction 
is given by ) ~ 7 ! 1 ^ 0 ) = + V nV'iV’-n)|0). When both states are occupied, the wavefunction is 

|^'l,2) ~7ll^l) =7l7ll^'o) = + V nV'iV'-n)|0)- As the hybridization between the bound states or the 

exchange coupling increases, the lower energy state becomes occupied and the ground state is I'I'i). When both states 
are below the chemical potential the ground state then becomes To determine the total energy of the system, 

one can diagonalize the Hamiltonian using a Bogoliubov transformation, H = ^n{(k) {ihln — where e„ is the 

energy of state n. The ground state energies are therefore®’^^ 

EA(^) = -\Y.\en{9)\. (Al) 


Appendix B: Weak Coupling Limit 


To obtain the energetically favorable magnetic orientation in the weak coupling limit, we solve Eq. (2) of the main 
text for the in-gap energies and expand to second order in l/kpr. We find the bound state energies are 


£ 71 ( 9 ) ~ An -1- Ao(A„ -|- Sn COS 9) 



(Bl) 


with 


Al 

Bl 

A2 

B2 


—2afa2(l — a|) -I- 2 aia 2 (l + af — 2afa2) cos2kFr 
{1 + + af) - ajil + aj)] 

—2a2ai{l + af — a2 — aia2){l — o;^) -I- 2a2ai{l — af + — afaf) cos2kFr 

(1 -k af)'^[af{l + af) - af(l -k af)] 

2 Q;f0:2(1 — af) — 20102(1 + af — 20^0!) cos 2kFr 
(1 -k al)'^[al(l + af) - af(l + af)] 

20102(1 — af + af — Oia2)(l — Oi) — 20102(1 -I- Oi — 02 — afaf) cos 2kF'i' 

{1 + af)'^[af{l + af) - af{l + af)] 


(B2) 


We consider three cases: when the bare energies are both above the chemical potential, both below the chemical 
potential, or on opposite sides of the chemical potential. The total energy of the system, according to Eq. (5) of the 
main text, is given by 


I -[ei(0) + e2(0)]/2, ei >0, £2 >0 

£i9) = { [e2{9)-ei{9)]/2, ei < 0, £2 > 0 . (B3) 

I [£i { 9 ) + £ 2 { 9)]/2 , £1 < 0, £2 < 0 


In all cases, the total energy is extremized when 0 = 0, tt, and for no intermediate values of 9. To determine the 
energetically favored magnetic configuration, we calculate 5£ = f(0) — £{'k). When both energies are above the 
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chemical potential, 


^ 

— = 2aia2 
^0 


( — 
ykpr 


1 + af + a2 + 2aia2 + afa^ + a^af + afa^ 

(l + a!y(l + aiy(l-a!al) 


1 — — a 2 ~ — afa^ — + afa^ 


= 2 g;iG ;2 


(l + af)^(l + ai)^(l-ofo|) 

(-L 

ykpr 


2 r 


1 I 4 4 
1 + 


il + aiy{l + aini-aial) 

ai+a2 + 2 a^a 2 + o;fa2 + 


cos2kpr 

(1 — cos2kpr) 


{l + air{l + aiy{l-alal) 

/ 1 X 2 

\kpr 


-(1 + cos2kpr) + 


4Q;fa2 


{l + aiy{l + aiy{l-alal) 


cos2kpr 


= 2cXiCX2 




[(1 + af)2(l + alY(l - alal) 
af + a2 + 4afQ;2 + afaf + afa^ 


(l + afy(l + al)^(l-afal) 
because ai, q ;2 < 1. Analogously, when ei, £2 < 0, 


(1 + cos2kpr) 


(1 — cos2kpr) 


> 0 


(B4) 


= —2aia2 
^0 


(— 

\kpr 


{l-afafy 


il + alni + afni-afaf) 


(1 — cos2kpr) 


a'i + af + Aafaf + afaf + d\a\ 

(l + a?) 2 (l + a 2 ) 2 (l_^ 2 ^ 2 ) 


(1 + cos2kpr) 


> 0 


(B5) 


because ai, 02 > 1 so that the preferred magnetic orientation is antiferromagnetic when the energies are on the same 
side of the chemical potential. Now suppose Ci, — 62 > 0, then we get 


S£ 


— 2 q;iQ ;2 




af + af + 2 ajQ ;2 + af + af + afaf + afaf 


(l + af)^(l + afy(af-af) 
af + af + 6afaf — af — af + afaf + afaf 


cos2kpr 


= 2aia2 


(l + af)^(l + af)^(af-af) 

/ 1 \ ^ af + af + 2 ajQ ;2 + afaf + afaf 
Iw [ (l + afy(l + afy(af-af) 

„,4 I „,4 

-(1 + cos2kpr) — 


(1 — cos2kpr) 


Aafaf 


(l + af)Hl + afnaf-af) 


= 2 q;iq ;2 


(af-af) 


af + a 2 + 4afaf + afaf + afaf 


(l + af)Hl + afr(af-af) 

.-2^4 


cos2kFr 


(l + afra + afnaf-af) 


(1 — cos2kpr) 


2\2 


{l + afyil + af)^{af-af) 


(1 + cos2kpr) 


< 0 


(B6) 


because 02 > ai- Therefore, making a similar argument when ei < 0 and £2 > 0, when the bare energies are on 
opposite sides of the chemical potential and sufficiently well separated, the impurities prefer to be oriented ferromag- 
netically. 

In the special case when ai = a 2 = a, the energy levels diverge according to Eq. (B2). The expansion of the bound 
state energies is instead given by 


,(d) « + (-l)"4a2Ao 

+ a^An 1 2a^ 


I cos(0/2)| sinkpr 
1 + a^ kpr 
1 — (1 — 2(4^) cos 2 kpr 

(l-a2)2(l+^2)3 


1 + q;'* — (1 — 4a^ + a^) cos 2 kpr 

(l-a 2 ) 2 (l + £,2)3 


cosO 


( — 
\kpr 


(B7) 


Although the leading order term contribution is of order exp{—r/^)/kpr and oscillates with 27r periodicity in 6 , the 
difference in total energy between the parallel and antiparallel conhgurations, when o; < 1 (cr > 1), again reduces to 
Eq. (B4) [Eq. (B5)] upon taking ai, a 2 —t a. 



Appendix C: Strong Hybridization Expressions 


For strongly hybridized identical impurities, the equations for the bound state energies are Eqs. (7) and (8), where 

2 


a = a < a 


1 + 


( — 
ykpr 


cos 2 kpr 


- 1 


b = 


siiikpr I 2 


c = a 


kpr 
2 + 


ykpr 


- 1 


(^) 


{cos2kpr — 1) 


d = a'^{l + 2 cos2kpr) 


( — 
\kpr 


(Cl) 


Appendix D: Additional Numerics 

We plot the difference in energies between ferromagnetic and antiferromagentic configurations 6 £ and S£gr for 
the same parameters as in Fig. 4 of the main text, i.e. a lattice of size x Ny = 33 x 25 with the chemical 
potential fi/t = 1, the superconducting gap Ag/t = 0.1, see Fig. 6. The exchange interaction strength is chosen to be 
J/t = 4. The difference between the YSR state energies in the collinear magnetic configurations is positive for nearly 
all values of r, indicating the antiferromagnetic orientation is preferred, which agrees with the analytics in the weak 
hybridization picture. The magnitude of the oscillations becomes larger if the quasiparticle contributions is included 
but the gap kept constant = Aq. Upon including renormalization of the gap, 6 £gr is significantly increased, similar 
to J/t = 1, again emphasizing the importance of the gap renormalization. 
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FIG. 6. The energy difference between ferromagnetic and 
antiferromagnetic configurations. The parameters are the 
same as in Fig. 4 of the main text with J/t = 4. 


FIG. 7. The sum of two YSR state energy 
found self-consistently. The parameters are 
the same as in Fig. 5 of the main text with 
J/f = 2.17. 


Taking J/t = 2.7 while leaving all other parameters the same as in Fig. 5 of the main text, we plot £ as a function 
of 0. Similar to the analytics, we find a ground state and a metastable state at the collinear configurations of the 
magnetizations. However, because of the self-consistent renormalization of the gap, there is a jump in £ at 0 « 7r/6 
where the QPT occurs. We note that the change in £ is several orders of magnitude smaller as compared with £gr, 
see Fig. 5 of the main text. 

















